library(Seurat)
library(cowplot)
library(umap)
library(dplyr)
library(Matrix)
library(readxl)
library(openxlsx)
library(ggplot2)
library(ggrepel)
library(ggpubr)
library(sctransform)
library(knitr)
library(kableExtra)
#library(biomaRt)
library(DESeq2)
library(escape)
library(dittoSeq)
library(GSEABase)
library(scater)
library(ComplexHeatmap)
gene.lists <- read_xlsx("Munoz_Yilmaz_CellCycle_signatures.xlsx")
gene.lists.names <- colnames(gene.lists)
GOI.lists <- c()
for (i in gene.lists.names) {
tmpList <- gene.lists %>% dplyr::select(all_of(i))
tmpList <- tmpList[!is.na(tmpList)]
GOI.lists[[i]] <- tmpList
}
the initial processing was done with r 3.6.1 with Seurat 3.2.0 - the UMAP comes out slightly differently in r 4.0.3 with Seurat 3.2.3*
#al.dat <- Read10X("200218Yil_data/al/filtered_feature_bc_matrix/")
#f.dat <- Read10X("200218Yil_data/f/filtered_feature_bc_matrix/")
#rf.dat <- Read10X("200218Yil_data/rf/filtered_feature_bc_matrix/")
#rf.rapa.dat <- Read10X("200218Yil_data/rf.rapa/filtered_feature_bc_matrix/")
#al <- CreateSeuratObject(counts = al.dat, project = "al", min.cells = 3, min.features = 200)
#f <- CreateSeuratObject(counts = f.dat, project = "f", min.cells = 3, min.features = 200)
#rf <- CreateSeuratObject(counts = rf.dat, project = "rf", min.cells = 3, min.features = 200)
#rf.rapa <- CreateSeuratObject(counts = rf.rapa.dat, project = "rf.rapa", min.cells = 3, min.features = 200)
#al[["percent.mt"]] <- PercentageFeatureSet(al, pattern = "^mt-")
#f[["percent.mt"]] <- PercentageFeatureSet(f, pattern = "^mt-")
#rf[["percent.mt"]] <- PercentageFeatureSet(rf, pattern = "^mt-")
#rf.rapa[["percent.mt"]] <- PercentageFeatureSet(rf.rapa, pattern = "^mt-")
#VlnPlot(al, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
#VlnPlot(f, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
#VlnPlot(rf, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
#VlnPlot(al, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
#merged_Seu <- merge(al, c(f,rf,rf.rapa), project = "Diet")
#merged_Seu <- NormalizeData(merged_Seu, normalization.method = "LogNormalize", scale.factor = 10000)
#merged_Seu <- FindVariableFeatures(merged_Seu, selection.method = "vst", nfeatures = 2000)
#merged_Seu <- ScaleData(merged_Seu)
#merged_Seu <- RunPCA(merged_Seu, features = VariableFeatures(object = merged_Seu))
#merged_Seu <- RunUMAP(merged_Seu, reduction="pca",dims=1:30)
#merged_Seu <- RunTSNE(merged_Seu, reduction="pca",dims=1:30)
#merged_Seu <- FindNeighbors(merged_Seu, dims = 1:30, verbose = FALSE)
#merged_Seu <- FindClusters(merged_Seu, verbose = FALSE)
#DimPlot(merged_Seu,reduction="umap",group.by="orig.ident",label=TRUE,repel=FALSE)
#FeaturePlot(merged_Seu,reduction="umap",features="mt-Cytb",min.cutoff=0,max.cutoff=4,cols=c("grey","red"))
#FeaturePlot(merged_Seu,reduction="umap",features="percent.mt",cols=c("grey","red"))
#FeaturePlot(merged_Seu,reduction="umap",features="nFeature_RNA",cols=c("grey","red"))
#save(merged_Seu, file="merged.Robj")
#load("merged.Robj")
#di <- SplitObject(merged_Seu, split.by = "orig.ident")
#for (i in 1:length(di)) {
# di[[i]] <- NormalizeData(di[[i]], verbose = FALSE)
# di[[i]] <- FindVariableFeatures(di[[i]], selection.method = "vst", nfeatures = 2000,
# verbose = FALSE)
#}
#dat.anchors <- FindIntegrationAnchors(object.list=di,dims=1:30)
#integrated <- IntegrateData(anchorset=dat.anchors,dim=1:30)
#DefaultAssay(integrated)<-"integrated"
#integrated <- ScaleData(integrated)
#integrated <- RunPCA(integrated,npcs=30)
#integrated <- RunUMAP(integrated,reduction="pca",dims=1:30)
#integrated <- RunTSNE(integrated,reduction="pca",dims=1:30)
#integrated <- FindNeighbors(integrated, dims = 1:30, verbose = FALSE)
#integrated <- FindClusters(integrated, verbose = FALSE)
NOTE: loading final object to avoid recalculating ssGSEA data
#load("integrated_orig.Robj")
load("../repo_data/integrated_final.Robj")
DefaultAssay(integrated)<-"integrated"
DimPlot(integrated,reduction="umap",split.by="orig.ident",group.by="orig.ident")
DimPlot(integrated,reduction="umap",group.by="orig.ident")
DimPlot(integrated,reduction="umap",group.by="integrated_snn_res.0.8", label=TRUE)
Using `as.character()` on a quosure is deprecated as of rlang 0.3.0.
Please use `as_label()` or `as_name()` instead.
integrated@meta.data$treat_clust <- paste0(integrated@meta.data$orig.ident,integrated@meta.data$integrated_snn_res.0.8)
integrated@meta.data$clust_treat <- paste0(integrated@meta.data$integrated_snn_res.0.8,integrated@meta.data$orig.ident)
integrated@meta.data$celltype <- integrated@meta.data$integrated_snn_res.0.8
Idents(object = integrated) <- integrated$celltype
integrated <- RenameIdents(object = integrated, '16' = '16_Tuft')
integrated <- RenameIdents(object = integrated, '11' = '11_EC')
integrated <- RenameIdents(object = integrated, '13' = '13_EEC')
integrated <- RenameIdents(object = integrated, '2' = '2_Stem')
integrated <- RenameIdents(object = integrated, '5' = '5_Stem')
integrated <- RenameIdents(object = integrated, '10' = '10_Stem')
integrated <- RenameIdents(object = integrated, '14' = '14_Paneth')
integrated <- RenameIdents(object = integrated, '8' = '8_Secretory_Progenitor')
integrated <- RenameIdents(object = integrated, '9' = '9_Goblet')
integrated <- RenameIdents(object = integrated, '1' = '1_Secretory_Progenitor')
integrated <- RenameIdents(object = integrated, '4' = '4_Secretory_Progenitor')
integrated <- RenameIdents(object = integrated, '15' = '15_Secretory_Progenitor')
integrated <- RenameIdents(object = integrated, '3' = '3_EC_Progenitor')
integrated <- RenameIdents(object = integrated, '6' = '6_EC_Progenitor')
integrated <- RenameIdents(object = integrated, '0' = '0_Early_TA')
integrated <- RenameIdents(object = integrated, '7' = '7_Early_TA')
integrated <- RenameIdents(object = integrated, '12' = '12_Unknown')
integrated[["clust_celltype"]] <- Idents(object = integrated)
Idents(object = integrated) <- integrated$clust_celltype
celltype_colors <- c("2_Stem"="#117733",
"5_Stem"="#999933",
"10_Stem"="#009E73",
"0_Early_TA"="#E69F00",
"7_Early_TA"="#D55E00",
"6_EC_Progenitor"="#0072B2",
"3_EC_Progenitor"="#56B4E9",
"11_EC"="#88CCEE",
"13_EEC"="#6699CC",
"1_Secretory_Progenitor"="#661100",
"4_Secretory_Progenitor"="#882255",
"8_Secretory_Progenitor"="#CC6677",
"15_Secretory_Progenitor"="#AA4499",
"14_Paneth"="#332288",
"16_Tuft"="#000000",
"9_Goblet"="#F0E442",
"12_Unknown"="#999999")
dp.cb <- DimPlot(integrated,reduction="umap", cols=celltype_colors, label=TRUE, repel=TRUE, pt.size=2, label.size=6) + NoLegend()
dp.cb
#pdf('Fig3b.cbSafe.pdf',width=14, height=10)
#dp.cb
#dev.off()
DefaultAssay(integrated)<-"RNA"
Idents(object = integrated) <- integrated$integrated_snn_res.0.8
plotOrder <- c("5","2","10","0","1","3","4","6","7","8","9","11","12","13","14","15","16")
vln_colors <- c("2"="#117733",
"5"="#999933",
"10"="#009E73",
"0"="#E69F00",
"1"="#661100",
"3"="#56B4E9",
"4"="#882255",
"6"="#0072B2",
"7"="#D55E00",
"8"="#CC6677",
"9"="#F0E442",
"11"="#88CCEE",
"12"="#999999",
"13"="#6699CC",
"14"="#332288",
"15"="#AA4499",
"16"="#000000")
Idents(integrated) <- factor(Idents(integrated), levels= plotOrder)
vl.cb <- VlnPlot(integrated, cols=vln_colors, idents = , features = c("Lgr5"), pt.size = 0.5, slot="data")
vl.cb
#pdf('Fig3c_cbSafe.pdf',width=14, height=8)
#vl.cb
#dev.off()
p.integrated <- table(integrated@meta.data$integrated_snn_res.0.8,integrated@meta.data$orig.ident)
round(prop.table(p.integrated,2),3)
al f rf rf.rapa
0 0.169 0.191 0.135 0.178
1 0.096 0.101 0.135 0.119
2 0.082 0.116 0.104 0.102
3 0.099 0.091 0.098 0.101
4 0.125 0.078 0.082 0.061
5 0.076 0.080 0.084 0.072
6 0.076 0.074 0.065 0.077
7 0.056 0.040 0.096 0.055
8 0.078 0.053 0.034 0.037
9 0.028 0.046 0.043 0.060
10 0.034 0.032 0.044 0.055
11 0.026 0.056 0.033 0.029
12 0.039 0.028 0.024 0.032
13 0.005 0.006 0.009 0.009
14 0.005 0.005 0.006 0.007
15 0.002 0.002 0.005 0.005
16 0.004 0.002 0.004 0.001
DefaultAssay(integrated)<- "integrated"
Idents(object = integrated) <- integrated$integrated_snn_res.0.8
dd <- subset(integrated, idents = c("2", "5", "10"), downsample=100)
topvst <- head(VariableFeatures(dd), 500)
mat <- as.matrix(dd@assays$integrated@scale.data) #as.matrix(subset_dd@assays$integrated@scale.data)
mat <- mat[topvst,]
genes <- c(GOI.lists$Biton_S1_ISC.I, GOI.lists$Biton_S1_ISC.II, GOI.lists$Biton_S1_ISC.III)
labels <- c(rep('Biton ISC I', length(GOI.lists$Biton_S1_ISC.I)),
rep('Biton ISC II', length(GOI.lists$Biton_S1_ISC.II)),
rep('Biton ISC III', length(GOI.lists$Biton_S1_ISC.III)))
idxs <- which(genes %in% rownames(mat))
genes <- genes[idxs]
labels <- labels[idxs]
mat <- mat[genes,]
ht <- Heatmap(mat, column_names_side = 'top', row_names_gp = gpar(fontsize = 9), column_names_gp = gpar(fontsize = 12),
column_title = '', name = 'scaled data', column_labels = rep('', ncol(mat)),
column_split = factor(as.character(dd$integrated_snn_res.0.8), levels = c('5', '2', '10')),
row_split = labels, row_order = genes, #column_order = sort(colnames(mat)),
cluster_column_slices = FALSE,
top_annotation = HeatmapAnnotation(cluster = as.character(dd$integrated_snn_res.0.8)))
draw(ht)
#pdf('Fig3D.pdf',width=12, height=10)
#draw(ht)
#dev.off()
data <- as.data.frame(read_excel('cluster_2_5_10stem_gsea.xlsx', sheet = 'Sheet2'))
data
data$FDR <- data$`FDR q-val` + 0.001
data$Comparison <- gsub('\\.noNA\\.NES','',data$Comparison)
comparisons <- c('f5_v_al5', 'rf5_v_al5', 'rf.rapa5_v_al5')
sub_data <- data[which(data$Comparison %in% comparisons),]
cl5 <- ggplot(data=sub_data, aes(y=NAME, x=factor(Comparison, levels = comparisons), size=-log10(FDR), color=NES)) +
geom_point() +
scale_color_gradient2(midpoint=0, low="blue", mid="white",
high="red", space ="Lab", limits=c(-3,3))+
scale_size_continuous(range=c(1,6))+
ggtitle('Cluster 5 GSEA') + theme_classic() +
theme(legend.position="right", axis.text.x = element_text(angle = 90)) + ylab('Gene Set') + xlab('Comparison')
cl5
#pdf('Fig5F.pdf',width=4, height=4)
#cl5
#dev.off()
comparisons <- c('f2_v_al2', 'rf2_v_al2', 'rf.rapa2_v_al2')
sub_data <- data[which(data$Comparison %in% comparisons),]
cl2 <- ggplot(data=sub_data, aes(y=NAME, x=factor(Comparison, levels = comparisons), size=-log10(FDR), color=NES)) +
geom_point() +
scale_color_gradient2(midpoint=0, low="blue", mid="white",
high="red", space ="Lab", limits=c(-3,3))+
scale_size_continuous(range=c(1,6))+
ggtitle('Cluster 2 GSEA') + theme_classic() +
theme(legend.position="right", axis.text.x = element_text(angle = 90)) + ylab('Gene Set') + xlab('Comparison')
cl2
#pdf('FigS3Da.pdf',width=4, height=4)
#cl2
#dev.off()
comparisons <- c('f10_v_al10', 'rf10_v_al10', 'rf.rapa10_v_al10')
sub_data <- data[which(data$Comparison %in% comparisons),]
cl10 <- ggplot(data=sub_data, aes(y=NAME, x=factor(Comparison, levels = comparisons), size=-log10(FDR), color=NES)) +
geom_point() +
scale_color_gradient2(midpoint=0, low="blue", mid="white",
high="red", space ="Lab", limits=c(-3,3))+
scale_size_continuous(range=c(1,6))+
ggtitle('Cluster 10 GSEA') + theme_classic() +
theme(legend.position="right", axis.text.x = element_text(angle = 90)) + ylab('Gene Set') + xlab('Comparison')
cl10
#pdf('FigS3Db.pdf',width=4, height=4)
#cl10
#dev.off()
#there is a glitch in this plot, cl2 loses x axis, legend order is different
figS3D <- ggarrange(cl2,cl10, ncol=2, nrow=1)
figS3D
#pdf('FigS3D_incorrect.pdf',width=8, height=4)
#figS3D
#dev.off()
DefaultAssay(integrated)<-"RNA"
colorScheme <- c("#C1BEBF","#fe6100")
fp.Muc2 <- FeaturePlot(integrated,reduction="umap",features="Muc2",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Muc2
fp.Tff3 <- FeaturePlot(integrated,reduction="umap",features="Tff3",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Tff3
fp.Lyz1 <- FeaturePlot(integrated,reduction="umap",features="Lyz1",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Lyz1
fp.Defa30 <- FeaturePlot(integrated,reduction="umap",features="Defa30",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Defa30
fp.Chga <- FeaturePlot(integrated,reduction="umap",features="Chga",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Chga
fp.Reg3a <- FeaturePlot(integrated,reduction="umap",features="Reg3a",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Reg3a
fp.Alpi <- FeaturePlot(integrated,reduction="umap",features="Alpi",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Alpi
fp.Atoh1 <- FeaturePlot(integrated,reduction="umap",features="Atoh1",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Atoh1
fp.Lgr5 <- FeaturePlot(integrated,reduction="umap",features="Lgr5",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Lgr5
fp.Smoc2 <- FeaturePlot(integrated,reduction="umap",features="Smoc2",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Smoc2
figS4A <- ggarrange(fp.Muc2,fp.Tff3,fp.Lyz1,fp.Defa30,fp.Chga,fp.Reg3a,fp.Alpi,fp.Atoh1,fp.Lgr5,fp.Smoc2, legend = "right", ncol = 5, nrow = 2)
figS4A
#pdf('FigS4A_cbSafe.pdf',width=20, height=10)
#figS4A
#dev.off()
DefaultAssay(integrated)<-"RNA"
gene.names <- rownames(integrated@assays$RNA@data)
Idents(object = integrated) <- integrated$seurat_clusters
stem.subset <- subset(integrated, idents = c("2","5","10"))
levels(stem.subset) <- c("5","2","10")
vln_colors <- c("2"="#117733",
"5"="#999933",
"10"="#009E73")
s.InData <- intersect(gene.names,GOI.lists$mm.s)
g2m.InData <- intersect(gene.names,GOI.lists$mm.g2m)
stem.subset[["percent.mm.s"]] <- PercentageFeatureSet(stem.subset, features = s.InData)
stem.subset[["percent.mm.g2m"]] <- PercentageFeatureSet(stem.subset, features = g2m.InData)
mm.s <- VlnPlot(stem.subset, cols = vln_colors, features="percent.mm.s",pt.size = 0.3,slot = "data")
mm.s
mm.g2m <- VlnPlot(stem.subset, cols = vln_colors, features="percent.mm.g2m",pt.size = 0.3,slot = "data")
mm.g2m
figS4b <- ggarrange(mm.s,mm.g2m, legend = FALSE, ncol=2, nrow=1)
figS4b
pdf('figS4b_cbSafe.pdf',width=8, height=4)
figS4b
dev.off()
png
2
integrated[["percent.mm.s"]] <- PercentageFeatureSet(integrated, features = s.InData)
integrated[["percent.mm.g2m"]] <- PercentageFeatureSet(integrated, features = g2m.InData)
mm.g2m.Flag <- if_else(integrated@meta.data$percent.mm.g2m >= 0.3, "Yes", "No")
integrated@meta.data$mm.g2m.Flag <- mm.g2m.Flag
p <- table(integrated$clust_treat,integrated$mm.g2m.Flag)
p.g2m.summary <- round(prop.table(p,2),3)
mm.s.Flag <- if_else(integrated@meta.data$percent.mm.s >= 0.2, "Yes", "No")
integrated@meta.data$mm.s.Flag <- mm.s.Flag
p <- table(integrated$clust_treat,integrated$mm.s.Flag)
p.s.summary <- round(prop.table(p,1),3)
This code is no longer working due to R and package updates but resulting data is stored in the seurat object escape 1.0.0 is probably required but is no longer available. the version in this R image is 1.0.1 escape ssGSEA - run ssGSEA to quantify expression of the BitonI gene set clusters
egs <- GeneSet(GOI.lists$Biton_S1_ISC.I, setName="BitonI")
ES <- enrichIt(obj = integrated, gene.sets = egs, groups = 1000, cores = 4)
[1] "Using sets of 1000 cells. Running 19 times."
Error in (function (classes, fdef, mtable) :
unable to find an inherited method for function ‘gsva’ for signature ‘"matrix", "character"’
# egs <- GeneSet(GOI.lists$Biton_S1_ISC.II, setName="BitonII")
# ES <- enrichIt(obj = integrated, gene.sets = egs, groups = 1000, cores = 4)
# integrated@meta.data$BitonII_ssGSEA <- ES$BitonII
# egs <- GeneSet(GOI.lists$Biton_S1_ISC.III, setName="BitonIII")
# ES <- enrichIt(obj = integrated, gene.sets = egs, groups = 1000, cores = 4)
# integrated@meta.data$BitonIII_ssGSEA <- ES$BitonIII
Figure 5G - Biton 1 in cluster5 Figure 5G_II - Biton II in cluster2 Figure 5G_III - Biton III in cluster10 Figure 5H - Pdgfa in cluster5 Figure S3E - Gkn3 in cluster5*
Idents(object = integrated) <- integrated$seurat_clusters
clust5.subset <- subset(integrated, idents = c("5"))
Idents(object = clust5.subset) <- clust5.subset$treat_clust
vln_colors <- c("al5"="#e69f00",
"f5"="#56b4f9",
"rf5"="#117733",
"rf.rapa5"="#d55e00")
plot.order <- c("al5","f5","rf5","rf.rapa5")
Idents(object = clust5.subset) <- factor(Idents(object = clust5.subset), levels = plot.order)
vln.BitonI.ssGSEA <- VlnPlot(clust5.subset, cols = vln_colors, features="BitonI_ssGSEA",slot = "data",pt.size = 0.4) +
stat_summary(fun = median, geom='point', size = 15, colour = "black", shape = 95) + NoLegend() + ylab("Enrichment Score") + ggtitle("Biton ISC-I_ssGSEA")
vln.BitonI.ssGSEA
#pdf('Fig5G.pdf',width=4, height=4)
#vln.BitonI.ssGSEA
#dev.off()
vln.Pdgfa <- VlnPlot(clust5.subset, cols = vln_colors, features = c("Pdgfa"), slot= "data",pt.size = 0.4) + NoLegend()
vln.Pdgfa
pdf('FigS4g_cbSafe.pdf',width=4, height=4)
vln.Pdgfa
dev.off()
png
2
vln.Gkn3 <- VlnPlot(clust5.subset, cols = vln_colors, features = c("Gkn3"), slot= "data",pt.size = 0.4) + NoLegend()
vln.Gkn3
pdf('FigS4h_cbSafe.pdf',width=4, height=4)
vln.Gkn3
dev.off()
png
2
clust2.subset <- subset(integrated, idents = c("2"))
Idents(object = clust2.subset) <- clust2.subset$treat_clust
vln_colors <- c("al2"="#e69f00",
"f2"="#56b4f9",
"rf2"="#117733",
"rf.rapa2"="#d55e00")
plot.order <- c("al2","f2","rf2","rf.rapa2")
Idents(object = clust2.subset) <- factor(Idents(object = clust2.subset), levels = plot.order)
vln.BitonII.ssGSEA <- VlnPlot(clust2.subset, cols = vln_colors, features="BitonII_ssGSEA",slot = "data",pt.size = 0.4) +
stat_summary(fun = median, geom='point', size = 15, colour = "black", shape = 95) + NoLegend() + ylab("Enrichment Score") + ggtitle("Biton ISC-II_ssGSEA")
vln.BitonII.ssGSEA
#pdf('Fig5G_II.pdf',width=4, height=4)
#vln.BitonII.ssGSEA
#dev.off()
clust10.subset <- subset(integrated, idents = c("10"))
Idents(object = clust10.subset) <- clust10.subset$treat_clust
vln_colors <- c("al10"="#e69f00",
"f10"="#56b4f9",
"rf10"="#117733",
"rf.rapa10"="#d55e00")
plot.order <- c("al10","f10","rf10","rf.rapa10")
Idents(object = clust10.subset) <- factor(Idents(object = clust10.subset), levels = plot.order)
vln.BitonIII.ssGSEA <- VlnPlot(clust10.subset, cols = vln_colors, features="BitonIII_ssGSEA",slot = "data",pt.size = 0.4) +
stat_summary(fun = median, geom='point', size = 15, colour = "black", shape = 95) + NoLegend() + ylab("Enrichment Score") + ggtitle("Biton ISC-III_ssGSEA")
vln.BitonIII.ssGSEA
#pdf('Fig5G_III.pdf',width=4, height=4)
#vln.BitonIII.ssGSEA
#dev.off()
bitonPlot <- ggarrange(vln.BitonI.ssGSEA,vln.BitonII.ssGSEA,vln.BitonIII.ssGSEA,nrow=1,ncol=3)
bitonPlot
#pdf('FigS4f_cbSafe.pdf',width=10, height=5)
#bitonPlot
#dev.off()
DefaultAssay(integrated)<-"RNA"
Idents(object = integrated) <- integrated$seurat_clusters
clust5.subset <- subset(integrated, idents = c("5"))
Idents(object = clust5.subset) <- clust5.subset$treat_clust
levels(clust5.subset) <- c("al5","f5","rf5","rf.rapa5")
clust5.subset$treat_clust <- factor(x = clust5.subset$treat_clust, levels = c("al5","f5","rf5","rf.rapa5"))
VlnPlot(clust5.subset, features = c("Oat"), split.by = "treat_clust", group.by = "treat_clust", slot= "data",pt.size = 0.5)
The default behaviour of split.by has changed.
Separate violin plots are now plotted side-by-side.
To restore the old behaviour of a single split violin,
set split.plot = TRUE.
This message will be shown once per session.
vln.Oat.5 <- VlnPlot(clust5.subset, features = c("Oat"), slot= "data",pt.size = 0.1)
vln.Oat.5
#pdf('Fig5H_Oat.5.pdf',width=4, height=4)
#vln.Oat.5
#dev.off()
clust2.subset <- subset(integrated, idents = c("2"))
Idents(object = clust2.subset) <- clust2.subset$treat_clust
levels(clust2.subset) <- c("al2","f2","rf2","rf.rapa2")
clust2.subset$treat_clust <- factor(x = clust2.subset$treat_clust, levels = c("al2","f2","rf2","rf.rapa2"))
VlnPlot(clust2.subset, features = c("Oat"), split.by = "treat_clust", group.by = "treat_clust", slot= "data",pt.size = 0.5)
vln.Oat.2 <- VlnPlot(clust2.subset, features = c("Oat"), slot= "data",pt.size = 0.1)
vln.Oat.2
#pdf('Fig5H_Oat.2.pdf',width=4, height=4)
#vln.Oat.2
#dev.off()
clust10.subset <- subset(integrated, idents = c("10"))
Idents(object = clust10.subset) <- clust10.subset$treat_clust
levels(clust10.subset) <- c("al10","f10","rf10","rf.rapa10")
clust10.subset$treat_clust <- factor(x = clust10.subset$treat_clust, levels = c("al10","f10","rf10","rf.rapa10"))
VlnPlot(clust10.subset, features = c("Oat"), split.by = "treat_clust", group.by = "treat_clust", slot= "data",pt.size = 0.5)
#vln.Oat.10 <- VlnPlot(clust10.subset, features = c("Oat"), slot= "data",pt.size = 0.1)
#vln.Oat.10
#pdf('Fig5H_Oat.10.pdf',width=4, height=4)
#vln.Oat.10
#dev.off()
DefaultAssay(integrated)<-"RNA"
Idents(object = integrated) <- integrated$seurat_clusters
clust5.2.10.subset <- subset(integrated, idents = c("5","2","10"))
Idents(object = clust5.2.10.subset) <- clust5.2.10.subset$treat_clust
levels(clust5.2.10.subset) <- c("al5","f5","rf5","rf.rapa5","al2","f2","rf2","rf.rapa2","al10","f10","rf10","rf.rapa10")
clust5.2.10.subset$treat_clust <- factor(x = clust5.2.10.subset$treat_clust, levels = c("al5","f5","rf5","rf.rapa5","al2","f2","rf2","rf.rapa2","al10","f10","rf10","rf.rapa10"))
VlnPlot(clust5.2.10.subset, features = c("Oat"), split.by = "treat_clust", group.by = "treat_clust", slot= "data",pt.size = 0.5)
vln.Oat.5.2.10 <- VlnPlot(clust5.2.10.subset, features = c("Oat"), slot= "data",pt.size = 0.1)
vln.Oat.5.2.10
#pdf('Fig5H_Oat.5.2.10.pdf',width=12, height=4)
#vln.Oat.5.2.10
#dev.off()
DefaultAssay(integrated)<-"RNA"
Idents(object = integrated) <- integrated$orig.ident
levels(integrated) <- c("al","f","rf","rf.rapa")
integrated$orig.ident <- factor(x = integrated$orig.ident, levels = c("al","f","rf","rf.rapa"))
VlnPlot(integrated, features = c("Oat"), split.by = "orig.ident", group.by = "orig.ident", slot= "data",pt.size = 0.5)
vln.Oat.noclust <- VlnPlot(integrated, features = c("Oat"), slot= "data",pt.size = 0.1)
vln.Oat.noclust
#pdf('Fig5H_Oat.noclust.pdf',width=4, height=4)
#vln.Oat.noclust
#dev.off()
# DefaultAssay(integrated)<-"RNA"
# Idents(object = integrated) <- integrated$treat_clust
# levels(x = integrated)
# ##
# #adjust the id1 and id2 variables to set up different tests
# #could use a loop for this
# ##
# id1 <- "al5"
# id2 <- "al2"
# outData <- paste(id1,id2,"Markers",sep=".")
# outData <- FindMarkers(integrated, ident.1 = id1, ident.2 = id2,logfc.threshold = 0, assay="RNA")
# outData <- tibble::rownames_to_column(outData,"#MGISym")
# rnk.tmp <- outData %>% dplyr::select(c("#MGISym","avg_logFC"))
# #write.table(rnk.tmp, sep='\t',file=paste0(id1,"_v_",id2,".rnk"),col.names=TRUE, quote=FALSE, row.names=FALSE)
#
# id1 <- "al10"
# id2 <- "al2"
# outData <- paste(id1,id2,"Markers",sep=".")
# outData <- FindMarkers(integrated, ident.1 = id1, ident.2 = id2,logfc.threshold = 0, assay="RNA")
# outData <- tibble::rownames_to_column(outData,"#MGISym")
# rnk.tmp <- outData %>% dplyr::select(c("#MGISym","avg_logFC"))
# #write.table(rnk.tmp, sep='\t',file=paste0(id1,"_v_",id2,".rnk"),col.names=TRUE, quote=FALSE, row.names=FALSE)
#
# id1 <- "al10"
# id2 <- "al5"
# outData <- paste(id1,id2,"Markers",sep=".")
# outData <- FindMarkers(integrated, ident.1 = id1, ident.2 = id2,logfc.threshold = 0, assay="RNA")
# outData <- tibble::rownames_to_column(outData,"#MGISym")
# rnk.tmp <- outData %>% dplyr::select(c("#MGISym","avg_logFC"))
# #write.table(rnk.tmp, sep='\t',file=paste0(id1,"_v_",id2,".rnk"),col.names=TRUE, quote=FALSE, row.names=FALSE)
#save(integrated, file="integrated_final.Robj")
sessionInfo()
R version 4.0.3 (2020-10-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.5 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 LC_PAPER=en_US.UTF-8
[8] LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] grid parallel stats4 stats graphics grDevices utils datasets methods base
other attached packages:
[1] ComplexHeatmap_2.6.2 scater_1.18.6 SingleCellExperiment_1.12.0 GSEABase_1.52.1 graph_1.68.0 annotate_1.68.0
[7] XML_3.99-0.5 AnnotationDbi_1.52.0 dittoSeq_1.2.6 escape_1.0.1 DESeq2_1.30.1 SummarizedExperiment_1.20.0
[13] Biobase_2.50.0 MatrixGenerics_1.2.1 matrixStats_0.58.0 GenomicRanges_1.42.0 GenomeInfoDb_1.26.7 IRanges_2.24.1
[19] S4Vectors_0.28.1 BiocGenerics_0.36.1 kableExtra_1.3.1 knitr_1.31 sctransform_0.3.2 ggpubr_0.4.0
[25] ggrepel_0.9.1 ggplot2_3.3.3 openxlsx_4.2.3 readxl_1.3.1 Matrix_1.2-18 dplyr_1.0.4
[31] umap_0.2.7.0 cowplot_1.1.1 Seurat_3.2.3
loaded via a namespace (and not attached):
[1] reticulate_1.18 tidyselect_1.1.0 RSQLite_2.2.3 htmlwidgets_1.5.3 BiocParallel_1.24.1 Rtsne_0.15 munsell_0.5.0
[8] codetools_0.2-16 ica_1.0-2 future_1.21.0 miniUI_0.1.1.1 withr_3.0.0 colorspace_2.0-0 rstudioapi_0.13
[15] ROCR_1.0-11 ggsignif_0.6.0 tensor_1.5 listenv_0.8.0 labeling_0.4.2 GenomeInfoDbData_1.2.4 polyclip_1.10-0
[22] farver_2.0.3 pheatmap_1.0.12 bit64_4.0.5 parallelly_1.23.0 vctrs_0.6.5 generics_0.1.0 xfun_0.21
[29] R6_2.5.1 clue_0.3-58 ggbeeswarm_0.6.0 rsvd_1.0.3 msigdbr_7.2.1 locfit_1.5-9.4 bitops_1.0-6
[36] spatstat.utils_2.0-0 cachem_1.0.3 DelayedArray_0.16.3 assertthat_0.2.1 promises_1.1.1 scales_1.1.1 beeswarm_0.2.3
[43] gtable_0.3.0 beachmat_2.6.4 Cairo_1.5-12.2 globals_0.14.0 goftest_1.2-2 rlang_1.1.4 genefilter_1.72.1
[50] GlobalOptions_0.1.2 splines_4.0.3 rstatix_0.6.0 lazyeval_0.2.2 broom_0.7.4 yaml_2.2.1 reshape2_1.4.4
[57] abind_1.4-5 backports_1.2.1 httpuv_1.5.5 tools_4.0.3 ellipsis_0.3.2 RColorBrewer_1.1-2 ggridges_0.5.3
[64] Rcpp_1.0.6 plyr_1.8.6 sparseMatrixStats_1.2.1 zlibbioc_1.36.0 purrr_1.0.2 RCurl_1.98-1.2 rpart_4.1-15
[71] openssl_2.2.0 deldir_0.2-9 GetoptLong_1.0.5 viridis_0.5.1 pbapply_1.4-3 zoo_1.8-8 haven_2.3.1
[78] cluster_2.1.0 magrittr_2.0.3 data.table_1.13.6 RSpectra_0.16-0 scattermore_0.7 circlize_0.4.12 lmtest_0.9-38
[85] RANN_2.6.1 fitdistrplus_1.1-3 GSVA_1.38.2 hms_1.0.0 patchwork_1.1.1 mime_0.9 evaluate_0.24.0
[92] xtable_1.8-4 rio_0.5.16 shape_1.4.5 gridExtra_2.3 compiler_4.0.3 tibble_3.0.6 KernSmooth_2.23-17
[99] crayon_1.4.1 htmltools_0.5.8.1 mgcv_1.8-33 later_1.1.0.1 tidyr_1.1.2 geneplotter_1.68.0 DBI_1.1.1
[106] MASS_7.3-53 car_3.0-10 cli_3.6.3 igraph_1.2.6 forcats_0.5.1 pkgconfig_2.0.3 foreign_0.8-80
[113] scuttle_1.0.4 plotly_4.9.3 xml2_1.3.2 vipor_0.4.5 webshot_0.5.2 XVector_0.30.0 rvest_0.3.6
[120] stringr_1.4.0 digest_0.6.36 RcppAnnoy_0.0.18 spatstat.data_2.0-0 rmarkdown_2.6 cellranger_1.1.0 leiden_0.3.7
[127] edgeR_3.32.1 uwot_0.1.10 DelayedMatrixStats_1.12.3 curl_5.2.1 shiny_1.6.0 rjson_0.2.20 lifecycle_1.0.4
[134] nlme_3.1-149 jsonlite_1.8.8 BiocNeighbors_1.8.2 carData_3.0-4 limma_3.46.0 viridisLite_0.3.0 askpass_1.1
[141] pillar_1.4.7 lattice_0.20-41 fastmap_1.2.0 httr_1.4.2 survival_3.2-7 glue_1.4.2 zip_2.1.1
[148] spatstat_1.64-1 png_0.1-7 bit_4.0.4 stringi_1.5.3 blob_1.2.1 BiocSingular_1.6.0 memoise_2.0.1
[155] irlba_2.3.3 future.apply_1.7.0
writeLines(capture.output(sessionInfo()), "2024_sessionInfo.txt")